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SUMMARY 

The global asymptotic nonlinear behavior of 11 explicit and implicit time discretizations for four 2x2 
systems of first-order autonomous nonlinear ordinary differential equations (ODEs) is analyzed. The 
objectives are to gain a basic understanding of the difference in the dynamics of numerics between the scalars 
and systems of nonlinear autonomous ODEs and to set a baseline global asymptotic solution behavior of 
these schemes for practical computations in computational fluid dynamics. We show how “numerical" basins 
of attraction can complement the bifurcation diagrams in gaining more detailed global asymptotic behav ior 
of time discretizations for nonlinear differential equations (DEs). We show how in the presence of spurious 
asymptotes the basins of the true stable steady states can be segmented by the basins of the spurious stable 
and unstable asymptotes. One major consequence of this phenomenon which is not commonly known is that 
this spurious behavior can result in a dramatic distortion and, in most cases, a dramatic shrinkage and 
segmentation of the basin of attraction of the true solution for finite time steps. Such distortion, shrinkage 
and segmentation of the numerical basins of attraction will occur regardless of the stability of the spurious 
asymptotes, and will occur for unconditionally stable implicit linear multistep methods. In other words, for 
the same (common) steady-state solution the associated basin of attraction of the DE might be very different 
from the discretized counterparts and the numerical basin of attraction can be very different from numerical 
method to numerical method. The results can be used as an explanation for possible causes of error, and slow r 
convergence and nonconvergence of steady-state numerical solutions when using the time-dependent 
approach for nonlinear hyperbolic or parabolic PDEs. 

KEY WORDS: Spurious steady-state numerical solutions, spurious asymptotes, global asymptotic 
behavior, nonlinear ODEs, numerical methods, time discretizations. 
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1. INTRODUCTION 

The tool that is utilized for the current study belongs to a multidisciplinary field of 
study in numerical analysis, sometimes referred to as “The Dynamics of Numerics 1 ”. 
Here the phrase “to study the dynamics of numerics” (dynamical behavior of a numeri- 
cal scheme) is restricted to the study of local and global asymptotic behavior and 
bifurcation phenomena of the nonlinear difference equations resulting from finite 
discretizations of a nonlinear differential equation (DE) subject to the variation of 
discretized parameters such as the time step, grid spacing, numerical dissipation 
coefficient, etc. In this paper, standard terminologies of nonlinear dynamics, chaotic 
dynamics (Guckenheimer and Holmes, 1983; Hale and Kocak, 1991) and computa- 
tional fluid dynamics (CFD) are assumed. For an introduction to the dynamics of 
numerics and its implications for algorithm development in CFD, see Yee et al. (1991) 
and Yee (1991) and references cited therein. 


1.1 Background 

The phenomenon that a nonlinear DE and its discretized counterpart can have 
different dynamical behavior (asymptotic behavior) was not uncovered fully until 
recently. Aside from truncation error and machine round-off error, a more fundamental 
distinction between the DE (continuum) and its discretized counterparts for genuinely 
nonlinear behavior is extra solutions in the form of spurious stable and unstable 
asymptotes that can be created by the numerical method. Here we use the term 
“discretized conterparts” to mean the finite difference equations (or discrete maps) 
resulting from finite discretizations of the underlying DE. Also we use the term “spurious 
asymptotic numerical solutions” to mean asymptotic solutions that satisfy the dis- 
cretized counterparts but do not satisfy the underlying ordinary differential equations 
(ODEs) or partial differential equations (PDEs). Asymptotic solutions here include 
steady-state solutions (fixed points of period one for the discretized equations), periodic 
solutions, limit cycles, chaos and strange attractors. See Section III and Guckenheimer 
and Holmes (1983), Hale and Kocak (1991) and Yee et al. (1991) for definitions. 

Iserles (1988) showed that while linear multistep methods (LM Ms) for solving ODEs 
possess only the fixed points (fixed points of period one) of the original DEs, popular 
Runge-Kutta methods may exhibit additional, spurious fixed points. It has been 
demonstrated by the authors and collaborators (Yee et al ., 1991; Yee, 1991; Sweby 
et al ., 1990; Griffiths et u/., 1992; Yee and Sweby, 1993a, 1993b) for nonlinear ODEs, 
and Lafon and Yee (1991, 1 992) for nonlinear reaction-convection model equations that 
such spurious fixed points as well as spurious fixed points of higher periods may be 
stable below the linearized stability limit of the scheme, depending on the initial data. 
Iserles et al. (1990), Hairer et al. (1989) and Humphries (1991) further advanced some 
theoretical understanding of the dynamics of numerics for initial value problems of 
ODEs. Iserles et al. and Hairer et al. classified and gave guidelines and theory on the 
types of Runge-Kutta methods that do not exhibit spurious period one or period two 
fixed points. Humphries (1991) showed that under appropriate assumptions if stable 


1 Named after the First IMA Conference on Dynamics of Numerics and Numerics of Dynamics, 
University of Bristol, England. July 31 August 2, 1990. 
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spurious fixed points exist as the time-step approaches zero, then they must either 
approach a true fixed point or become unbounded. However, convergence in practical 
calculations involves a finite time step At as the number of integrations n oo rather 
than Af->0, as oo. There appear to be missing links between theoretical develop- 
ment and practical scientific computation. Our aim is to provide some of these missing 
links that were not addressed in Iserles (1988), Iserles et al (1990), Hairer et al (1989), 
Humphries (1991) and our earlier work. In particular, we want to show in more detail 
the global asymptotic behavior of time discretizations when finite but not extremely 
small Ar is used. Other aspects that were not addressed in Iserles (1988) for different 
iteration procedures in solving the resulting nonlinear algebraic equations are reported 
in greater depth in our companion papers (Yee and Sweby, 1993a, 1993b). 

L2 Relevance and Motivations 

Although the understanding of the dynamics of numerics of systems of nonlinear ODEs 
and PDEs is important in its own right and has applications in the various nonlinear 
scientific fields, our main emphasis is CFD applications. Time-marching types of 
methods (time-dependent approach) are commonly used in CFD because the steady 
PDEs of higher than one dimension are usually of the mixed type. When a time- 
dependent approach is used to obtain steady-state numerical solutions of a fluid flow 
or a steady PDE, a boundary value problem (BVP) is transformed into an initial- 
boundary value problem (IBVP) with unknown initial data. If the steady PDE is 
strongly nonlinear and/or contains stiff nonlinear source terms, phenomena such as 
slow convergence, nonconvergence or spurious steady-state numerical solutions and 
limit cycles commonly occur even though the time step is well below the linearized 
stability limit and the initial data are physically relevant. One of our goals is to search 
for logical explanations for these phenomena via the study of the dynamics of numerics. 
Here the term “time-dependent approach" is used loosely to include some of the 
iteration procedures (due to implicit time discretizations), relaxation procedures, and 
preconditioners for convergence acceleration strategies used to numerically solve 
steady PDEs. This is due to the fact that most of these procedures can be viewed as 
approximations of time-dependent PDEs (but not necessarily the original PDE that 
was under consideration). If one is not careful, numerical solutions other than the 
desired one of the underlying PDE can be obtained (in addition to spurious asymptotes 
due to the numerics). 

One consequence of the existence of stable and unstable spurious asymptotes below 
or above the linearized stability limit of the numerical schemes is that these spurious 
features may greatly affect the dynamical behavior of the numerical solution in practice 
due to the use of a finite time step. As discussed in details in later sections and also in 
Yee et al. (1991), Yee and Sweby (1993a, b), Lafon and Yee (1992), Sweby and Yee 
(1991), Yee et al (1992), it is possible that for the same steady-state solution, the 
associated basin of attraction of the underlying DEs (which initial conditions lead to 
which asymptotic states) might be very different from that of the basin of attraction of 
the discretized counterparts due to the existence of spurious stable and unstable 
asymptotic numerical solutions. In other words, there is a separate dependence on 
initial data for the individual DEs and their discretized counterparts. Here the basin of 
attraction is a domain of a set of initial conditions whose solution curves (trajectories) 
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all approach the same asymptotic state. Also we use the term “exact” and “numerical” 
basins of attraction to distinguish “basins of attraction of the underlying DEs” and 
“basins of attraction of the discretized counterparts”. 

In view of the spurious dynamics, it is possible that numerical computations may 
converge to an incorrect steady state or other asymptote which appears to be physically 
reasonable. One major implication is that what is expected to be physical initial data 
associated with the underlying steady state of the DE might lead to a wrong steady 
state, a spurious asymptote, or a divergence or nonconvergence of the numerical 
solution. In addition, the existence of spurious limit cycles may result in the type of 
nonconvergence of steady-state numerical solutions observed in time-dependent 
approaches to the steady states. It is our belief that the understanding of the symbiotic 
relationship between the strong dependence on initial data and permissibility of 
spurious stable and unstable asymptotic numerical solutions at the fundamental level 
can guide the tuning of the numerical parameters and the proper and/or efficient usage 
of numerical algorithms in a more systematic fashion. It can also explain why certain 
schemes behave nonlinearly in one way but not another. Here strong dependence on 
initial data means that for a finite time step At that is not sufficiently small, the 
asymptotic numerical solutions and the associated numerical basins of attraction 
depend continuously on the initial data. Unlike nonlinear problems, the associated 
numerical basins of attraction of linear problems are independent of At as long as At is 
below a certain upper bound. 

Nonunique Steady-State Solutions of Nonlinear DEs vs. Spurious Asymptotes : The phe- 
nomenon of generating spurious steady-state numerical solutions (or other spurious 
asymptotes) by certain numerical schemes is often confused with the nonuniqueness (or 
multiple steady states) of the DE. In fact, the existence of nonunique steady-state 
solutions of the continuum can complicate the numerics tremendously (e.g., the basins 
of attraction) and is independent of the occurrence of spurious asymptotes of the 
associated scheme. But, of course, a solid background in the theory of nonlinear ODEs 
and PDEs and their dynamical behavior is a prerequisite in the study of the dynamics of 
numerics for nonlinear PDEs. See Yee et a/., 1991 for a discussion. It is noted that the 
approach and primary goal of our work is quite different from the work of e.g., Beam 
and Bailey (1988) and Jameson (1991). The main goal of Beam and Bailey (1988) and 
Jameson (1991) was to study the nonunique steady-state solutions admitted by the PDE 
as the physical parameter is varied. Our primary interest is to establish some working 
tools and guidelines to help delineate the true physics from numerical artifacts via the 
dynamics of numerics approach. The knowledge gained from our series of studies (Yee 
et al. , 1991; Lafonand Yee, 1991; Lafon and Yee, 1992) hopefully can shed some light on 
the controversy about the existence of multiple steady-state solutions through numeri- 
cal experiments for certain flow types of the Euler and/or Navier Stokes equations. 

1.3 Objectives and Outline 

The primary goal of the series of papers (present and the companion papers Yee et al . , 
1991; Yee and Sweby, 1993a, b; Lafon and Yee, 1991; Lafon and Yee, 1992) is to lay the 
foundation for the utilization of the dynamics of numerics in algorithm development for 
computational sciences in general and CFD in particular. This is part II of this series of 
papers on the same topic. Part I (Yee et al. , 1991) concentrated on the dynamical behavior 
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of time discretizations for scalar nonlinear ODEs. The intent of part I was to serve 
as an introduction to motivate this concept to researchers in the field of CFD and to 
present new results for the dynamics of numerics for first-order scalar autonomous 
ODEs. 

The present paper, the second of this series, is devoted to the study of the dynamics 
of numerics for 2 x 2 systems of ODEs. Here we show how “numerical” basins of 
attraction can complement the bifurcation diagrams in gaining more detailed global 
asymptotic behavior of numerical methods for nonlinear DEs. We show how in the 
presence of spurious asymptotes the basins of the true stable steady states can be 
segmented by the basin of the spurious stable and unstable asymptotes. One major 
consequence of this phenomenon which is not commonly known is that this spurious 
behavior can result in a dramatic distortion and, in most cases, a dramatic shrinkage 
and segmentation of the basin of attraction of the true solution for finite time steps. 
Such distortion, shrinkage and segmentation of the numerical basins of attraction 
will occur regardless of the stability of the spurious asymptotes, and will occur for 
unconditionally stable implicit linear multistep methods. In other words, for the same 
steady-state solution, the associated basin of the DE might be very different from its 
discretized counterparts. The basins can also be very different from numerical method 
to numerical method. The present study reveals for the first time the detail interlock- 
ing relationship of numerical basins of attraction and the causes of error, and slow 
convergence and nonconvergence of steady-state numerical solutions when using the 
time-dependent approach. 

The article of Lafon and Yee ( 1 99 1 ), the third of this series, was devoted to the study 
of the dynamics of numerics of commonly used numerical schemes in CFD for a model 
reaction-convection equation. The article of Lafon and Yee (1992), the fourth of 
this series, was devoted to a more detailed study of the effect of numerical treatment 
of nonlinear source terms on nonlinear stability of steady-state numerical solution 
for the same model nonlinear reaction-convection BVP. In our companion papers 
(Sweby et ai , 1990; Griffiths et al, 1991a, 1992b), a theoretical bifurcation analysis of 
a class of explicit Runge-Kutta methods and spurious discrete travelling wave phenom- 
enon were presented. In yet another companion paper, Yee and Sweby (1993a), the global 
asymptotic nonlinear behavior of three standard iterative procedures in solving 
nonlinear systems of algebraic equations arising from four implicit LMMs is analyzed 
numerically. 

I A Outline 

The outline of this paper is as follows. Section II discusses the connection of the 
dynamics of numerics for systems of ODEs and numerical approximations of time- 
dependent PDEs. Section III reviews background material for nonlinear ODEs and 
their numerical methods. Section IV describes four 2x2 systems of nonlinear 
first-order autonomous model ODEs. Section V describes the 1 1 time discretizations 
and the associated bifurcation diagrams for the four model ODEs. Section VI discusses 
the combined basins of attraction and bifurcation diagrams for the underlying 
schemes. Comparison between a linearized implicit Euler and Newton method is 
briefly discussed in Section 6.5. The paper ends with some concluding remarks in 
Section VII. 
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2. THE DYNAMICS OF NUMERICS OF SYSTEMS OF ODEs AND 
NUMERICAL APPROXIMATIONS OF TIME-DEPENDENT PDEs 

For finite discretizations of PDEs, spurious asymptotes and especially spatially- 
varying spurious steady states can be independently introduced by time and spatial 
discretizations (Yee et al , 1991; Lafon and Yee, 1991; Lafon and Yee, 1992). The interaction 
between temporal and spatial dynamical behavior is more complicated when one is 
dealing with the nonseparable temporal and spatial finite-difference discretizations 
such as the Lax-Wendroff type. The analysis and the different features of the numerics 
due to temporal and spatial discretizations can become more apparent by separable 
temporal and spatial finite difference methods (FDM). A standard method for obtain- 
ing such a FDM is the method of lines (MOL) procedure where the time-dependent 
PDE is reduced to a system of ODEs (by replacing the spatial derivatives by finite differ- 
ence approximations). The resulting approximation is called semi-discrete, since the 
time variable is left continuous. The semi-discrete system in turn can be solved by the 
desired time discretizations. Similar semi-discrete systems can be obtained by finite 
element methods except in this case an additional mass matrix is involved. Besides the 
MOL approach, coupled nonlinear ODEs can arise in many other ways when analyzing 
nonlinear PDEs. See for example Globus et al (1991), Hung et al (1991), Foias et al 
(1985), Temam (1989), Kwak (1991), Schecter and Shearer (1990), and Shearer et al 
(1987). Among these possibilities, the idea of inertial manifold (IM) and approximate 
inertial manifold (AIM) for incompressible Navier-Stokes (Foias et al 1985; Temam, 
1989; Kwak, 1991), the relationship between shock waves, heteroclinic orbits of systems 
of ODEs (Schecter and Shearer, 1990; Shearer et al. , 1987), and flow visualization of 
numerical data (Globus et al , 1991; Hung, 1991) are touched upon here. 

2.1 Asymptotic Analysis of the Method of Lines Approach 

When the ODEs are obtained from a semi-discrete approximations of PDEs, the 
resulting system of ODEs contains additional system parameters and discretized 
parameters as opposed to physical problems governed by ODEs. Depending on the 
number of grid points “J” used, the dimensions of the resulting system of semi-discrete 
approximations of ODEs can be very large. Also, depending on the differencing scheme 
the resulting discretized counterparts of a PDE can be nonlinear in A /, the grid spacing 
Ax and the numerical dissipation parameters, even though the DEs consist of only one 
parameter or none. One major consideration is that one might be able to choose 
a “safe” numerical method to solve the resulting system of ODEs to avoid spurious 
stable steady states due to time discretizations. However, spurious steady states and 
especially spatially varying steady states introduced by spatial discretizations in 
nonlinear hyperbolic and parabolic PDEs for CFD applications appear to be more 
difficult to avoid. In the case of the MOL approach, if spurious steady states due to 
spatial discretizations exist, the resulting ODE system has already inherited this 
spurious feature as part of the exact solutions of the semi-discrete case. We remark that 
spurious stable and unstable asymptotes other than the steady states due to time 
discretizations are also more difficult to avoid than spurious steady states. See Sections 
V and VI for some illustrations. Taking for example the nonlinear ODE models that are 
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considered, it is relatively easy to avoid spurious steady states due to time discretiz- 
ations since, if a numerical steady state U* for the ODEdU/dt = S(U) is spurious, then 
S(U*) / 0. This is not the case for spurious asymptotes such as limit cycles. 

In addition to the aforementioned considerations, it is well known from the theory 
of nonlinear dynamics for ODEs that much of the established theory and known 
behavior of nonlinear dynamics are restricted to lower dimensional first-order ODEs 
(or for problems that exhibit lower dimensional dynamical behavior). Moreover, if 
higher than two-time level numerical methods are used, the dynamics of these 
discretized counterparts usually are richer in structure and more complicated to 
analyze than their two-time level cousins. Therefore, in order to gain a first hand 
understanding of the subject we restrict our study to 2 x 2 systems of first-order 
autonomous ODEs and two-time level numerical methods with a fixed time step, even 
though the current study is far removed from the realistic setting. Studies of 3 x 3 
systems and general J x J systems are in progress. 

Due to the complexity of the subject matter, this paper concerns fixed time step (and 
fixed grid spacing) time-marching methods only. The fixed or local variable time step 
control method study can also shed some light on identifying whether certain flow 
patterns are steady or unsteady. See Yee et al (1990) for some examples. Proper 
regulation of a variable time step to prevent the occurrence of spurious steady-state 
numerical solutions will be a subject of future research. In order to isolate the different 
causes and cures of slow convergence and nonconvergence of time-marching methods, 
our study concerns nonlinearity and stiffness that are introduced by DEs containing 
smooth solutions. Nonlinearity and stiffness that are introduced by the scheme, the 
coupling effect in the presence of a source term (terms) in coupled system of PDEs, the 
highly stretched nonuniform structured and unstructured grids, the discontinuities in 
grid interfaces and/or the discontinuities inherent in the solutions, and external flows 
that need special boundary conditon treatment with a truncated finite computation 
domain are added factors and require additional treatment or different analysis. These 
are not considered at the moment. Generalization of our study to include grid adaption 
as one of the sources of nonlinearity and/or stiffness introduced by the numerics is 
reported in Sweby and Yee (1964) and Budd et al (1994). 

2.2 Inertial Manifold (IM ) and Approximate Inertial Manifold ( AIM ) 

The concept of IMs was introduced by Foias et al (1985). See Foias et al (1985), 
Temam (1989) and Kwak (1991) for details of the subject. The key idea of IMs and 
AlMs is to establish theories to aid in better understanding of nonlinear phenomena 
and turbulence via the study of the interaction of short and long wavelengths of 
dissipative systems. Basically, an I M is a finite-dimensional submanifold that contains 
all the attractors and invariant sets of an infinite-dimensional dynamical system 
described by some dissipative PDEs. It establishes the criterion for the reduction of 
long-term dynamics of certain infinite-dimensional problems to a finite system of 
ODEs. An attractive feature is that the reduction introduces no error in the problem. 
That is, the IM contains all pertinent information about the long-term dynamics of the 
original system. One of the main objectives of AIMs is to handle cases where the IMs 
are not known to exist. AIMs also can help in finding good algorithms for dealing with 
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the IMs that are known to exist. AIMs may also help reduce finite but extremely large 
systems of ODEs to lower-dimensional problems. In a nut shell, the derivation of IMs 
and AIMs is based on the decomposition of the unknown function into large-scale and 
small scale components. In the case of fluid dynamics, those structures can be identified 
as large and small eddies. Thus an IM or AIM corresponds to an exact or approximate 
interaction law between the short and long wavelengths. Kwak (1991 ) showed that the 
long-term dynamics of some two-dimensional incompressible Navier-Stokes equations 
can be completely described by a finite system of ODEs. Kwak does so by finding 
a nonlinear change of variable that embeds the incompressible Navier-Stokes equa- 
tions in a system of reaction-diffusion equations that possess an IM. All of the theories 
of IMs and AIMs are very involved and interested readers are encouraged to read Foias 
et ai (1985), Temam (1989) and Kwak (1991) and the references cited therein. 

2.3 Relationship Between Shock Waves and. Heteroclinic Orbits of Systems of ODEs 

Another example of the importance of understanding the “dynamics” and the “dy- 
namics of numerics” of systems of ODEs is related to the study of shocks using 
equilibrium bifurcation diagrams of associated vector fields. This was introduced by 
Shearer et ai (1987). The authors find of great interest how one can reduce the study of 
admissible shock wave solutions of a 2 x 2 hyperbolic conservation laws to the study of 
heteroclinic orbits of a system of nonlinear ODEs. Further development in this area 
can help in constructing suitable approximate Riemann solvers in numerical computa- 
tions. Schecter and Shearer (1990) studied undercompressive shocks for nonstrictly 
hyperbolic conservation laws by adding information to the equilibrium bifurcation 
diagrams (introduced by Shearer et ai) about heteroclinic orbits of the vector fields. 
The augmented equilibrium bifurcation diagrams are then used in the construction of 
solutions of Riemann problems. 

2 A Dynamics of Numerics and Flow Visualizations of Numerical Data 

The use of flow visualization of numerical data (numerical solutions of finite discretiz- 
ations of e.g., fluid flow problems) in an attempt to understand the true flow physics has 
become increasingly popular in the last decade. See, Globus et ai (1991) and Hungef ai 
(1991) and references cited therein. Many of the techniques rely on the extraction of 
the boundary surfaces by analyzing a set of appropriate vector fields. Approximations 
are then performed based on this set of vector fields. The study of the topological 
features of certain flow physics based on the numerical data is then related to the 
study of fixed points of the associated systems of ODEs. Fluid problems with known 
flow physics can be used to reveal how well the associated vector fields of the 
numerical data can mimic the true physics. It can also help to delineate spurious 
flow patterns that are solely due to the numerics. At the present time we are entering 
into the regime where CFD is extensively used to aid the understanding of complicated 
flow physics that is not amenable to analysis otherwise. In the situation where the 
numerical data indicate flow structures which are not easily understood, a good 
understanding of the spurious dynamics that can be introduced by the numerics is 
needed. 
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3. PRELIMINARIES 

Consider a 2 x 2 system of first-order autonomous nonlinear ODEs of the form 

dU 

j7 = S(1/), (3.1) 

where U and S are vector functions of dimension 2, and S(U) is nonlinear in U. 
A fixed point U E of an autonomous system (3.1) is a constant solution of (3.1); 
that is 


S(U E ) = 0, (3.2) 

where the subscript “E” stands for “exact" and U E denotes the fixed points of the ODE 
as opposed to the additional fixed points of the discretized counterparts (spurious fixed 
points) due to the numerical methods which we will encounter later. 

Let the eigenvalues of J{U E ) = (dS/dU)\ VE (the Jacobian matrix of S(U) evaluated 
at U E ) be and X 2 . Here J(U E ) is assumed to be nonzero. The fixed point U E is 
hyperbolic if Re(A t ) ^ 0, i = 1, 2. If both X k are real, U E is a saddle if A l X 2 < 0 and a node if 
X 1 X 2 >0. If exactly one 2 f = 0, then U E is semihyperbolic. If the eigenvalues are 
complex, then U E is a spiral. The “tightness” of the spiral is governed by the magnitude 
of the imaginary part of the eigenvalues. If the eigenvalues both have a zero real part, 
then U E is non-hyperbolic. Such a fixed point is called a center. Under this situation, 
more analysis is needed to uncover the real behavior of (3.1) around a non-hyperbolic 
fixed point. The fixed point U E is stable if both and X 2 have negative real parts. U E is 
unstable if a has a positive real part. In the non-hyperbolic case the fixed point is 
neutral. 

If due to a variation of a parameter of the ODE a fixed point becomes unstable, then, 
if at the point of instability the eigenvalues are distinct and real, the resulting bifurcation 
will be to another fixed point. Such bifurcation is called a steady bifurcation. If, however, 
the eigenvalues are complex, then the bifurcation will be of a Hopf type. This is a 
slightly simplified classification, since our main concern in this work is not on the 
variation of the ODE parameter. Detailed background information can be found in 
(Guckenheimer and Holmes, 1983; Hale and Kocak, 1991). 

Consider a nonlinear discrete map from a finite discretization of (3.1) 

U"+ i = U n + D(U\r\ (3.3) 

where r = At and D(U n , r) is linear or nonlinear in r depending on the numerical 
method. A fixed point U D of (3.3) is defined by £/" + 1 = or 


U D =U D + D(U D ,r) (3.4) 

or D(U Dl r) — 0. A fixed point V D of period p > 0 of (3.3) is defined by L/ n + p = U n 
with U n + k ^U n of k < p. In the context of discrete systems, the term “fixed point” 
without indicating the period means “fixed point of period 1” or the steady-state 
solution of (3.3). Here we use the term asymptote to mean a fixed point of any period, 
a limit cycle (in the discrete sense - invariant set), chaos, or a strange attractor. 
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The type of finite discretization of (3.1) represented in (3.3) assumed the use of 
two-time level schemes. Otherwise the vector dimension of (3.3) would be 2{k — \) 
instead of 2 where k is the number of the time level of the scheme. Here the vector 
function D is assumed to be consistent with the ODE (3.1) in the sense that fixed points 
of the ODE are fixed points of the scheme; however, the reverse need not hold. It is this 
feature accompanied by other added dynamics, that the discretized counterparts of the 
underlying ODE possess a much richer dynamical behavior than the original ODE. 
Thus the fixed points U D of D( U D , r) = 0 may be true fixed points U E of (3. 1 ) or spurious 
fixed points U s . The spurious fixed points U s are not roots of S(U) = 0. That is 
S(U s )^0. Spurious asymptotes are asymptotic numerical solutions of (3.3) but not 
(3.1). 

Letting U 9 = U D + S n , then a perturbation analysis on (3.3) yields 

,3.5, 

Assuming dO(U DJ r)/dU ^0, then the fixed point U D is stable if the eigenvalues of 
J D = / 4- dD(U Dy r)/dU lie inside the unit circle. If both eigenvalues are real and both lie 
inside (outside) the unit circle, then the fixed point is a stable (unstable) node. If one is 
inside the unit circle and the other outside, then the fixed point is a saddle. If both 
eigenvalues are complex, then the fixed point is a spiral. If the eigenvalues lie on the unit 
circle, then the fixed point of (3.3) is indeterminant and additional analysis is required to 
determine the true behavior of (3.3) around this type of fixed point. For a more refined 
definition and the difference in fixed point definition between ODEs and discrete maps, 
see Panov et al (1956), Perron (1929) and Hsu (1987) and references cited therein. The 
reader is referred to Guckenheimer and Holmes (1983), Hale and Kocak (1991), 
Langford and looss (1980), and Werner (1980) for full details on the subject of 
bifurcation theory. 

An important feature which can arise (for both systems of ODEs (3.1) and their 
discretizations) as the result of a Hopf bifurcation is a limit cycle where the trajectory 
traverses a closed curve in phase space. In all but a few simple cases such limit cycles are 
beyond analysis. 


4. MODEL 2x2 SYSTEMS OF NONLINEAR FIRST-ORDER 
AUTONOMOUS ODEs 

Four 2x2 systems of nonlinear first-order autonomous model ODEs are considered. 
The systems considered with U T = (u, r) or z = u -f iv are a 


1. Dissipative complex model: 


dz . 

— = z(i + r. - \z\ 2 ) 


(4.1) 
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2. Damped Pendulum model: 



du 

dr 


— ev — sin(u) 


3. Predator-Prey model: 


du 

dr 


3w + 4 u 2 — 0.5uu — u 3 


du 

— = — 2Av + uv 
dr 


(4.2a) 

(4.2b) 

(4.3a) 

(4.3b) 


4. Perturbed Hamiltonian System model: 
du 

^- = £(1 -3u)+|[l — 2u + u 2 — 2u(l -u)] (4.4a) 


— = £(1 — 3u) — | [1 —2v + v 2 — 2v(l — u)] (4.4b) 

Here e is the system parameter for (4.1), (4.2) and (4.4). 

The perturbed Hamiltonian model can be related to the numerical solution of the 
viscous Burgers’ equation with no source term 


du I d(u 2 ) 
St + 2~dx~ 



(4.5) 


Let Uj(t) represent an approximation to u(x p t) of (4.5) where Xj=jAxJ = 
with Ax the uniform grid spacing. Consider the three-point central difference in space 

assume Xj= i w, = constant, which implies that 
i dn ; /df — 0* If we take 7 = 3 and Ax = 1/3 then, with e = 9/?, this system can be 
reduced to a 2 x 2 system of first-order nonlinear autonomous ODEs (4.4) with 
U = (w 1? u 2 ) = (w, v). In this case, the nonlinear convection term is contributing to the 
nonlinearity of the ODE system (4.4). 


These four equations were selected to bring out the dynamics of numerics for four 
different types of solution behavior of the ODEs. The dissipative complex system (4.1) 
possesses either a unique stable fixed point or limit cycle with an unstable fixed point 
depending on the value of c. This is the rare situation where the analytical expression of 
a limit cycle can be found. The purpose of choosing (4.1) is to illustrate the numerical 


accuracy of computing a limit cycle and the spurious dynamics associated with this 
type of asymptote. The damped pendulum (4.2), arising from modelling of a physical 
process, exhibits a periodic structure of an infinite number of fixed points. The 
predator-prey model (4.3), arising from modelling of biological process, exhibits 
multiple stable fixed points without a periodic pattern as model (4.2). The perturbed 
Hamiltonian model (4.4), which arises as a gross simplification of finite discretization of 
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the viscous Burgers' equation, exhibits an unique stable fixed point. Following the 
classification of fixed points of (3.1) in Section III, one can easily obtain the following: 

Fixed Point of (4.1 ): The dissipative complex model has a unique fixed point at 
(u,e) = (0,0) for e ^0. The fixed point is a stable spiral if £ < 0. It is a center if £ = 0. For 
e > 0, the fixed point (0,0) becomes unstable with the birth of a stable limit cycle with 
radius equal to y/c centered at (0, 0). Figure 4.1 shows the phase portrait {u~v plane) of 
system (4. 1 ) for t: = — 1 and £ = 1 respectively. Here the entire (u, t?) plane belongs to the 
basins of attraction of the stable fixed point (0,0) if £ < 0. On the other hand, if £ > 0, the 
entire (u, v) plane except the unstable fixed point (0< 0) belongs to the basin of attraction 
of the stable limit cycle centered at (0,0). 

Fixed Points of (4.2): The damped pendulum (4.2) has an infinite number of fixed 
points, namely (kn, 0) for integer k. If k is odd, the eigenvalues of the Jacobian J(U E ) are 
of opposite sign and these fixed points are saddles. If k is even, however, two cases must 
be considered, depending on the value of c. If e < 2 and positive, the eigenvalues are 
complex with negative real part and the fixed points are stable spirals. If £ ^ 2, the 
eigenvalues are real and negative and the fixed points are nodes. If £ =0, the spirals 
become centers. Figure 4.2 shows the phase portrait and their corresponding basins of 
attraction for system (4.2). The different shades of grey regions represent the various 
basins of attraction of the respective stable fixed points for c = 0.5 and c = 2.5. 

Fixed Points of (4.3 ): The fixed points of the predator-prey equation are less regular 
than those for the damped pendulum equation. System (4.3) has four fixed points (0, 0), 
(0, 1), (3,0) and (2.1, 1.98). By looking at the eigenvalues of the Jacobian of S, one finds 
that (0,0) is a stable node, (2.1, 1.98) is a stable spiral, and (TO) and (3,0) are saddles. 
Figure 4.3 shows the phase portrait and their corresponding basins of attraction for 
system (4.3). The different shades of grey regions represent the various basins of 
attraction of the respective stable fixed points. The white region represents the basin of 




u u 


Figure 4.1 Phase Portraits and basins of Attraction Dissipative Complex Equation. 



DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 231 
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Figure 4.2 Phase Portraits and Basins of Attraction Damped Pendulum Equation. 




Figure 4.3 Phase Portraits and Basins of Attraction Predator- Prey Equation. 
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1.0 =3 0 = 3 




u u 

Figure 4.4 Phase Portraits and Basins of Attraction Viscous Burger’s Equation (Central Difference in 
Space). 


divergent solutions. Note that the trajectories near the unstable separatrices actually 
do not merge with the unstable branch of separatrices, but only appear to merge due to 
the thick drawings of the solution trajectories. 

Fixed Points of ( 4.4 ): The perturbed Hamiltonian has four steady-state solutions of 
which three are saddles and one is a stable spiral at (1/3, 1/3) for e^O. For e = 0 the 
stable spiral becomes a center. Figure 4.4 shows the phase portrait and their corre- 
sponding basins of attraction for system (4.4). The shaded region represents the basins 
of attraction for the fixed point (1/3, 1/3) for a = 0 and £ = 0.01. The white region 
represents the basin of divergent solutions. From here on we refer to (4.4) also as 
a viscous Burgers' equation with central difference in space. 


5. NUMERICAL METHODS AND BIFURCATION DIAGRAMS 

This section describes the 1 1 time discretizations and their corresponding bifurcation 
diagrams for the four model ODEs (4.1)— (4.4). The 1 1 numerical methods are listed in 
Section 5.1. Section 5.2 discusses the stability of selected fixed points of the discretized 
counterparts of the model ODEs as functions of system parameters. Section 5.3 
discusses the bifurcation diagrams as a function of the discretized parameter At with the 
system parameter held fixed. 

5.1 Numerical Methods 

The 9 explicit and two implicit methods considered are the explicit Euler, two 
second-order Runge-Kutta, namely, the modified Euler (R-K 2) and the improved 
Euler (R-K 2), two third-order Runge-Kutta (R-K 3), a fourth-order Runge-Kutta 
(R-K 4), the two and three-step predictor-corrector (Lambert, 1973), and noniterative 
linearized forms of the implicit Euler and the trapezoidal methods. 
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(1) Explicit Euler (lst-order; R-K 1): 

U" +i = U n + rS n - S" = S(U n ), (5.1) 

(2) Modified Euler (R-K 2): 

U" +i = U n + rs(u n + ^S n y ( 5 . 2 ) 

(3) Improved Euler (R-K 2): 

t/- +1 = (7" + ^[S" + S(l/ + rS")], (5.3) 

(4) Heun (R-K 3): 


W +l = U”+ r -(k 1+ 3k 3 ) (5.4) 

k 2 = S n 

k 2 = s(u" + r -k 
L k 3 = s(u n + jk 2 y 

(5) Kutta (R-K 3): 


U n+l = U n + r -(k l +4k 1 +k i ) 

k j =S n 

k 2 = + 

k 3 = S(V-rk l +2rk 2 ), 

(6) R-K 4: 


(5.5) 


V +l = V + ~(k l +2k 2 + 2k 3 + k 4 ) 
k t =S" 

k2 = s(u"+ r -k l ) 


(5.6) 
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/ k 3 = s(u" + r -k 2 ^ 

fc 4 = S( U" + rk 3 ), 

(7, 8) Predictor-corrector for m = 2, 3 (PC2, PC3): 


U t0) = U" + rS " 


U (k + 1 »=(/" + - [S" + S<“], /c = 0, 1, . . . , m - 1 


(yn+l = (;n+ _[ S n + 5 . m -l , ] ) 


(9) Adam-Bashforth (2nd-order): 


U" + ' = U n + - [3 S(t/") - S(U"~ *)], 


(5.7) 


(5.8) 


(10) Linearized Implicit Euler: 

L/" + 1 = U" + r{I-rJ n ) 'S" 


r = 



and det(/ - rJ n ) ¥= 0, 


(11) Linearized Trapezoidal: 


fjn+l 


U" + r | 


/- 



S” 


(5.9) 


(5.10) 


where the numeric identifier after the t4 R-K” indicates the order of accuracy of the 
scheme and r = A? and det( ) means the determinant of the quantity inside the ( ). 
Schemes (10) and (11) are unconditionally stable methods. See Beam and Warming 
(1976) and Yee (1989) for the versatility of the linearized implicit Euler and linearized 
trapezoidal methods in CFD applications. A comparison between Newton method in 
solving the steady part of the ODEs and the linearized implicit method (5.9) for model 
(4.4) is included in Section 6.5. Studies on Newton method in solving the steady state 
part of the PDE and some iteration procedures in solving the nonlinear algebraic 
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equation resulting from four implicit LMMs are reported in a separate paper (Yee and 
Sweby 1933a). Although the explicit Euler can be considered as an R-K 1, it is also 
a LMM. All of the R-K methods (higher than first order) and the predictor-corrector 
methods are nonlinear in the parameter space r, and all LMMs are linear in r. As 
discussed in Yee et al. (1991), a necessary condition for a scheme to produce spurious 
fixed points of period one is the introduction of nonlinearity in the parameter space r. It 
can be shown later that this property plays a major role on the shapes and sizes of the 
associated numerical basins of attraction of the scheme. For simplicity in referencing, 
hereafter we use “implicit Euler” and “trapezoidal” to mean the linearized forms (5.9) 
and (5.10), respectively, unless otherwise stated. 


5.2 Stability of Fixed Points oj Numerical Methods as a Function of System 
Parameters 

In our later study, we assume a fixed system parameter so that only the discretized 
parameter comes into play. However, in order to get a feel for the numerical stability of 
these schemes around selected stable fixed points U F as a function of the system param- 
eter e, Figures 5. 1-5. 3 show the stability regions of the schemes as a function of the 





Figure 5.1 Stability Regions vs. System Parameters Dissipative Complex Equation. 
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Figure 5,1 ( Continued ) 


system parameter e around a selected fixed point for each of the models. The linearized 
stability regions for the R-K methods of the same order behave in exactly the same 
manner, and the linearized stability regions around stable U E of the linearized 
implicit methods are not interesting, since they have the same regions of stability as 
the ODEs. 

The stability diagrams presented were obtained by numerically solving the absolute 
stability polynomials for the various methods, in most cases using Newton iteration. 
For the Runge-Kutta schemes (of order p ^ 4) the stability (Griffiths et al. y 1992; 
Lambert, 1973) condition is that 

2 p r p 

1 + Ar + — \ — <1, (5dl) 

p! 

where A are the eigenvalues of the Jacobian of S(l/). For the Predictor-Corrector of 
steps p = 2, 3 the stability condition is that 

:p r p 

1 + Ar + - ■ ■ + — — <1, 

2 P 


(5.11) 
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and for the Adams-Bashforth method the roots fj. of 



(5.13) 


satisfy | ^ | < 1 . Note that all of these expressions only hold for the U E fixed points of the 
system. 

In all cases the boundary of the stability region is when unit modulus is attained. The 
linearized implicit Euler and trapezoidal methods are unconditionally stable for the 
stable exact fixed points U E of the ODE systems we are considering. 

These stability regions can be used to isolate the key regions of the e parameter to be 
considered for the study of dynamics of numerics later. Due to the enormous number of 
possibilities, detailed study can only concentrate on one to two representative system 
parameters. Even with such a restriction, as can be seen later, computing the corre- 
sponding bifurcation diagrams and basins of attraction is very CPU intensive. Fortu- 
nately the computation can be made highly parallel. Figures 5. 1-5.3 also can serve as 
a spot check on the numerical results presented in the next section. 



Figure 5.2 Stability Regions vs. System Parameters Damped Pendulum Equation. 
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Figure 5.2 ( Continued ) 


5.3 Bifurcation Diagrams 

In this section, we show the bifurcation diagrams of selected R-K methods. It illus- 
trates some of the many ways in which the dynamics of a numerical discretization 
of 2 x 2 first-order autonomous nonlinear system of ODEs can differ from the 
system itself. Note that there is no limit cycle or higher dimensional tori counter- 
part for the scalar first-order autonomous ODEs. Spurious limit cycles and higher 
dimensional tori can only be introduced by the numerics when solving nonlinear 
ODEs other than scalar first-order autonomous ODEs (if 2-time level schemes are 
used) and/or by using a scheme with higher than two-time level for the scalar first- 
order autonomous ODEs. In Section VI, we showed how numerical basins of attrac- 
tion can complement the bifurcation diagrams in gaining more detailed global 
asymptotic behavior of numerical schemes. We purposely present our results in 
this order (not showing the basins of attraction) in order to bring out the importance 
of basins of attraction for the time-dependent approach in obtaining steady-state 
numerical solutions. 

Even though the analytical solutions of the these models are known, depending on 
the scheme, the dynamics of their discretized counterparts might be very difficult to 
analyze. In particular, some analytical linearized analysis (without numerical computa- 
tions) of fixed points of periods one and two is possible for the predator-prey and the 
damped pendulum case. However, analytical analysis for the dissipative complex 
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model and the perturbed Hamiltonian is not practical. For a detailed analysis of these 
selected cases, readers are referred to Sweby and Yee ( 199 1 ). For the majority of the cases 
where rigorous analysis is impractical we study the dynamics of numerics using numerical 
experiments. 

Note that some global solution behavior of fixed points of the nonlinear discretized 
equations (5.1)-(5.10) for (4.1)-(4.4) can be obtained by the pseudo arclength continua- 
tion method devised by Keller (1977), a standard numerical method for obtaining 
bifurcation curves in bifurcation analysis. A major shortcoming of the pseudo arclength 
continuation method is that for problems with complicated bifurcation patterns, it 
cannot provide the complete bifurcation diagram without known start up solutions 
for each of the main bifurcation branches before one can continue the solution along 
a specific main branch. For spurious asymptotes it is usually not easy to locate even just 
one solution on each of these branches. 

The nature of our calculations requires thousands of iterations of the same equation 
with different ranges of initial data on a preselected (u, v) domain and range of the 
discretized parameter space At. Since the NASA Ames CM-2 allows vast numbers 
(typically 65,536) of calculations to be performed in parallel, our problem is perfect for 
computation on the CM-2. 






Figure 5.3 Stability Regions vs. System Parameters Viscous Burgers' Equation (Central Difference in 
Space). 
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PC2 


PC3 




t 


Figure 5.3 ( Continued ) 


To obain a “full” bifurcation diagram, the domain of initial data and the range of the 
Af parameter are typically divided into 512 equal increments. For each initial datum 
and Af, the discretized equations are preiterated 3,000 5,000 (more or less depending 
on the ODE and scheme) before the next 4,000-6,000 iterations are plotted. The 
preiterations are necessary in order for the trajectories to settle to their asymptotic 
value. The high number of iterations plotted (overlay on the same plot) is to detect 
periodic orbits or invariant sets. Since the results are a three dimensional graph 
((Af, u, t 1 )), we have taken slices in a given constant v- and u-plane in order to enhance 
viewing the decrease CPU computations. Note that with this method of computing the 
bifurcation diagrams, only the stable branches are plotted. Some of the bifurcation 
diagrams in a v = constant plane for the four model ODEs and for the modified Euler, 
improved Euler, Kutta and R-K 4 methods are shown in Figures 5.4-5. 8. Figure 5.4 
shows a typical example of spurious stable fixed points (branches 3 and 4 on the 
diagram) occurring below the linearized stability by the modified Euler method. It also 
shows the existence of spurious asymptotes such as limit cycles, higher order periodic 
solutions and possibly numerical chaos (chaos introduced by numerics), See later 
sections and subsections for further details. Selected bifurcation diagrams for the rest of 
the numerical methods are illustrated in Section IV with basins of attraction 
superimposed (see Figures 6.3 6.5,6.13, 6.14, and 6.19 6.20). See also the original 
NASA internal report RNR-92-008, March 1992 for additional illustrations. Due to the 
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Un 



Figure 5.4 Bifurcation Diagram Predator-Prey Equation Modified Euler. 


plotting package, the labels (u„, v„) on all of the figures are («", v"). In the plots, r = At 
unless stated. 

The term full bifurcation as defined in Yee et al. (1991) is used to mean bifurcation 
diagrams that cover the essential lower-order periods in such a way as to closely resolve 
the true bifurcation diagram of the underlying discrete map for a selected range of 
initial data domains. This is necessary since solutions with different initial conditions 
will converge to different asymptotic limits. All of the computations shown are “full” 
bifurcation diagrams. 

The following summarizes the spurious dynamical behavior of the 1 1 numerical 
methods based on selected domains of initial data and ranges of the discretized 
parameter r. Numerical results agree with the analytical linearized analysis reported in 
Sw-eby and Yee (1991). 


Bifurcation Diagrams of Numerical Methods for Model (4.1 ): For /: = 0, (4.1) is nondis- 
sipative (or a Hamiltonian system), and all of the 1 1 numerical methods which are 
non-simplectic converge quite slowly to the fixed point (0,0). We conjecture that 
simplectic schemes (Sanz-Serna, 1990) would be more appropriate for /; = 0. For 
sufficiently small negative (positive) e, all of the studied schemes converge extremely 
slowly to the stable spiral (limit cycle). This is a typical example of slow convergence 
of the numerical solution due to the stiffness of the system parameter. While the 
bifurcation diagrams for *;^0 for the various numerical methods are not too 
interesting, the bifurcation diagrams for t; > 0 are very instructive. Figure 5.5 shows the 
bifurcation diagrams for the four R-K methods for r, = 1. 
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Figure 5.5 Bifurcation Diagrams Dissipative Complex Equation, c= 1, v = 0.0. 


Note also that R-K 4 method gives the most overall accurate numerical approxima- 
tions of the true limit cycle with radius ^/e centered at (0,0). The Adam-Bashforth, 
PC2, PC3, implicit Euler and trapezoidal methods give the least accurate numerical 
approximation of the limit cycle for r closer to the linearized stability. The R-K 4 and 
Heun methods produced spurious higher-order limit cycles (invariant set of multiple 
circles on the diagrams). See Section IV and Figures 6.6 and 6.8 for more details. These 
diagrams illustrate the unreliability of trying to compute a true limit cycle with any 
sizable r. This should not be surprising since the scheme only gives an 0(r p ) approxima- 
tion to the solution trajectories. In addition, since the limit cycle is not a fixed point, we 
would expect inaccuracies to be introduced. However, inaccuracies are not easy to 
detect in practice, especially when a numerical solution produces the qualitative 
features expected. See Section VI and Figures 6.3-6.9 for additional details. All of the 
studied explicit methods produce spurious asymptotes. 

For e > 0, the trapezoidal method produces no spurious steady states. However, the 
implicit Euler method in addition to maintaining an unconditionally stable feature of 
the exact limit cycle, also turns the unstable fixed point U E — (0, 0) of the ODE (4. 1 ) into 
a stable fixed point for r ^ 1. See Figures 6.5, 6.8 and 6.9 for additional details. 
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u n Modified Euler u n Improved Euler 
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Figure 5.6 Bifurcation Diagrams Damped Pendulum Equation, e = 1 , v = 0.0. 


Bifurcation Diagrams of Numerical Methods for Model (4.2): All of the studied explicit 
and implicit methods produce spurious asymptotes. In particular, some of the explicit 
methods (even explicit Euler) produce spurious limit cycles for certain e values. For 
certain ranges of r and ^ values the implicit Euler and trapezoidal methods turn the 
saddle points of (4.2) into an unstable fixed point of different type (see Figure 6. 1 2). For 
the modified Euler method, spurious steady states occur below the linearized stability 
limit of the scheme. See Section VI and Figures 6.10-6.12 for additional details. 

Bifurcation Diagrams of Numerical Methods for Model (4.3): Again, all of the 
studied explicit and implicit methods generate spurious asymptotes. Also, some of 
the explicit methods produce spurious limit cycles. For certain ranges of the r, the 
trapezoidal method turns the saddle points (exact fixed points of (4.3)) into unstable 
fixed points of different types, and the implicit Euler method turns the saddle points 
into stable fixed points of different type. The numerical results coincide with analytical 
analysis by examining the eigenvalues of the Jacobian of the resulting discrete map. 
Transcritical bifurcations introduced by the R-K 4 method resulted in the production 
of spurious steady state below (and very near) the linearized stability limit of the 
scheme. See, Section VI and Figures 6.13 6.18 for additional details. 



244 


H. C YEE AND P K SWEBY 


u n 



4 1 improved Euler 

3 —^ 

: 1 . 

-i i 1 1 1 1 1 1 1 1 1 r = aDt 

0 2 0.6 1.0 1.4 1.8 


u n 


Kutta 


3 



H 


0 



u n 



Figure 5.7 Bifurcation Diagrams Predator- Prey Equation, v — 0.0. 


More than one spurious fixed point below the linearized stability of the scheme was 
introduced by the modified Euler method (see Fig. 5.4). From the form of the modified 
Euler scheme it is easily seen that as well as the exact fixed points U E of the ODEs, any 
other value U s satisfying 


V s + ^S(U s )=U E (5.14) 

will also be a fixed point of the scheme. As mentioned earlier, we refer to these 
additional fixed points as spurious fixed points. Note that the U E on the right-hand side 
of (5. 14) encompasses both stable and unstable fixed points of the ODE and so, for the 
predator-prey equations (since S contains cubic terms in U), there are up to twelve (real) 
spurious steady states, three for each exact fixed point U E . In fact there are six such 
spurious steady states which lie in the v = 0 plane. All of them occur below the 
linearized stability limits of the exact fixed points, although not all are stable there. 
Four (stable ones) of the six are shown in the bifurcation diagram of Figure 5.4, 
numbered 1, 3, 4, 6 of the bifurcation branch. The other two are unstable. Note also that 
the branch numbered 6 is in fact not stable but represents the stable eigen-direction 
(separatrix) in the r = 0 plane of a saddle point. 
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Figure 5.8 Bifurcation Diagrams Viscous Burgers’ Equation. i: = 0.l, r = 0.333 (Central Difference in 
Space). 


Bifurcation Diagrams of Numerical Met hods for Model (4.4 ): For t: = 0, the ODE (4.4) is 
nondissipative and thus for small r, slow convergence was experienced. For r beyond 
the linearized limit and with e = 0 all of the explicit methods produce spurious limit 
cycles. For z > 0 (and not too large) all of the studied 1 1 explicit and implicit methods 
produce spurious asymptotes. Also, all of the explicit methods produce spurious limit 
cycles. For £ = 0.1, the Kutta and Heun methods introduce spurious asymptotes 
(higher than period one) that are below the linearized stability limit of the scheme. See 
Figures 6.19 6.25 for additional details. 


6 BASINS OF ATTRACTION AND BIFURCATION DIAGRAMS 

This section illustrates how basins of attraction can complement the bifurcation 
diagrams in gaining more detailed global asymptotic behavior of time discretizations 
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for nonlinear DEs. Two different representations of the numerical basins of attrac- 
tion were computed on the NASA Ames CM-2. One representation is bifurcation 
diagrams as a function of At with numerical basins of attraction superimposed on 
a constant r- or u-plane. The other representation is the numerical basins of attraction 
with stable asymptotes superimposed on the phase plane (u, i ). Before discussing 
numerical results for each of the model ODEs, the next subsection gives some 
preliminaries on how to compute and on how to interpret the basins of attraction 
diagrams for the CM-2. 

6.1 Introduction 

To obtain a bifurcation diagram with numerical basins of attraction superimposed on 
the CM-2, the preselected domain of initial data on a constant v- or u-plane and the 
preselected range of the At parameter are divided into 5 1 2 equal increments. Again, for 
the bifurcation part of the computations, with each initial datum and At, the discretized 
equations are preiterated 3,000 5,000 steps before the next 5,000 iterations (more or 
less depending on the problem and scheme) are plotted. The bifurcation curves appear 
on the figures as white curve, white dot and white dense dots. While computing the 
bifurcation diagrams it is possible to overlay basins of attraction for each value of At 
used. For the numerical basins of attraction part of the computation with each value of 
At used, we keep track of where each initial datum asymptotically approaches and 
color code them (appearing as a vertical strip) according to the individual asymptotes. 
While efforts were made to match color coding of adjacent strips on the bifurcation 
diagram, it was not always practical or possible. Care must therefore be taken when 
interpreting these overlays. 

For the basins of attraction on the phase plane (u,r) with selected values of At 
and the stable asymptotes superimposed, the (u, r) domain is divided into 512 x 512 
points of initial datum. With each initial datum and the selected At, we preiterate 
the respective discretized equation 3,000-5,000 steps and plot the next 5,000 steps 
to produce the asymptotes. Again, for the basins of attraction part of the compu- 
tations, for each value of At used, we keep track of where each initial datum 
asymptotically approaches and color code them according to the individual asym- 
ptotes. All of the selected time steps At shown are based on the bifurcation diagram 
with the basins of attraction superimposed. The chosen time steps were selected 
to illustrate special features of the different bifurcation phenomena on the (u,t*) 
plane. Details of the techniques used for detection of asymptotes and basins of 
attraction are given in the appendix of Sweby and Yee (1991). Note that in all of 
the plots, if color printing is not available, the different shades of grey represent different 
colors. 

As a prelimiary, and before discussing our major results, we discuss the numerical 
basins of attraction associated with modified Euler, improved Euler, Kutta and R-K 4 
methods for the two scalar first-order autonomous nonlinear ODEs studied in part I of 
our companion paper (Yee et at ., 1991). The two scalar ODEs are: 




( 6 . 1 ) 
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and 

du 

— = au(\ — u)(0.5 — u). (6.2) 

The fixed points for (6.1) with a > 0 are u = 0 (unstable) and u= 1 (stable), and no 
additional higher order periodic fixed points or asymptotes exist. The basin of attraction 
for the stable fixed point u = 1 is the entire positive plane for all values of a > 0. 

The fixed points for (6.2) with a > 0 are u = 0 (unstable), u = 1 (unstable) and u = 0.5 
(stable) and no additional higher-order periodic solutions or asymptotes exist. The 
basin of attraction for the stable fixed point u = 0.5 is 0 < u < 1 for all a > 0. The white 
curve, white dots and white dense dots of Figures 6.1 and 6.2 show the bifurcation 
diagrams for four of the R-K methods for (6.1) and (6.2). For more details of the dynamics 
of numerics for systems (6. 1 ) and (6.2), see Yee et al. (1991). Intuitively, in the presence of 
spurious asymptotes the basins of the true stable steady states can be separated by the 
numerical basins of attraction of the stable and unstable spurious asymptotes. 

Take, for example, the ODE (6. 1 ) where the entire domain u is divided into two basins 
of attraction for the ODE independent of any real a. Now if one numerically integrates 
the ODE, depending on the scheme and r, extra stable and unstable fixed points of any 
order can be introduced by the scheme. The bifurcation part of Figures 6.1 and 6.2, 
cannot distinguish the types of bifurcation and the periodicity of the spurious fixed 
points of any order. With the numerical basins of attraction and their respective 
bifurcation diagrams superimposed on the same plot, the type of bifurcation and to 
which initial data asymptote to which stable asymptotes become apparent. Note that 
for Figures 6. 1 and 6.2 r ~ aA t. 

For example, any initial data residing in the green region in Figure 6.1 for the 
modified Euler method belong to the numerical basin of attraction of the spurious 
(stable) branch emanating from u = 3 and r = 1 . Thus, if the initial data is inside the 
green region, the solution can never converge to the exact steady state using even 
a small fixed but finite Af (all below the linearized stability limit of the scheme). Note 
that the green region extends upward as r decreases below 1 . Thus for certain ranges of 
r values, the domain is divided into four basins (instead of two for the ODE). But of 
course higher period spurious fixed points exist for other ranges of r and more basins 
are created within the same u domain. 

A similar situation exists for the R-K 4 method (Fig. 6.1), except now the numerical 
basins of attraction of the spurious fixed points occur very near the linearized stability 
limit of the scheme, with a small portion occurring below the linearized stability limit. 
In constrast to the improved Euler method (Fig. 6.1), the green region represents the 
numerical basins of one of the spurious stable transcritical bifurcation branches of the 
fixed point. The bifurcation curve directly below it with the corresponding red portion 
is the basin of the other spurious branch. See Yeect al. (1991) or Hale and Kocak (1991) 
for a discussion of the different types of bifurcations. With this way of color coding the 
basins of attraction, one can readily see (from the plots) that for ODE (6. 1 ), the modified 
Euler, improved Euler and R-K 4 methods, experience one steady bifurcation before 
a period doubling bifurcation occurs (Fig. (6.1)). Using the PC3 method to solve (6.1) 
(figure not shown; see Yee et al ., 1991), more than one consecutive steady bifurcation 
occurs before period doubling bifurcation. For ODE (6.2), the improved Euler experi- 
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Figure 6.1 (See Color Plate 1 at the back of this issue ! 


ences two consecutive steady bifurcations before a period doubling bifurcation occurs 
(Fig. (6.2)). Using the PC3 method to solve (6.2)(figure not shown; see Yee cl ai 1991), 
four consecutive steady bifurcations occur before period doubling bifurcations. The 
modified Euler and R-K 4 methods, however, experience only one steady bifurcation 
before period doubling bifurcations occur. 
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Figure 6.2 (See Color Plate M at the back of this issue.) 


The next section presents similar diagrams for the 2 x 2 systems of model nonlinear 
ODEs (4.1) (4.4). In this case, only basins of attraction with bifurcation diagrams 
superimposed on v = constant planes are shown. Selected results for both representa- 
tions of numerical basins of attraction are shown in Figures 6.3 6.5, 6.8 6.25 for the 
numerical methods. Section 6.6 summarizes similar results presented in Yee and 
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Sweby (1993a) for iterative procedures in solving nonlinear systems of algebraic 
equations arising from four implicit LMMs. In the plots r = At. White dots and white 
curves on the basins of attraction with bifurcation diagrams superimposed represent 
the bifurcation curves. White dots and white closed curves on the basins of attraction 
with the numerical asymptotes superimposed represent the stable fixed points, stable 
periodic solutions or stable limit cycles. The black regions represent divergent 
solutions. 

Note that the streaks on some of plots are either due to the non-settling of the 
solutions within the prescribed number of preiterations or the existence of small 
isolated spurious asymptotes. Due to the high cost of computation, no further 
attempts were made to refine their detailed behavior since our purpose was to show 
how, in general, the different numerical methods behave in the context of nonlinear 
dynamics. 

6.2 Numerical results for the Dissipative Complex Equation 

Figures 6.3-6.S, 6.8, 6.9 show selected results for the two representations of numerical 
basins of attraction for model (4.1) for e = 1. The exact solution for (4.1) with e = 1 is 
a stable limit cycle with unit radius centered at (0, 0). The basin of attraction for the limit 
cycle is the entire (u,t;) plane except the unstable fixed point (0,0). 

Comparing Figures 6.3 6. 5 with Figure 5.5, one can appreciate the added informa- 
tion that the basin of attraction diagrams can provide. As At moves closer to the 
linearized stability limit of the limit cycle, the size (red) of the numerical basins of 
attraction decreases rapidly. This is due to the existence of spurious unstable asym- 
ptotes below as well as above the linearized stability limit. The green region, shown in 
Figure 6.5 using the implicit Euler method, is the numerical basin of attraction for the 
stabilized fixed point (0,0). Note how the implicit Euler method turns an unstable fixed 
point (0, 0) of the ODE system into a stable one for At ^ 1 . 

Figures 6.6 and 6.7 show the phase trajectories and Figures 6.8 and 6.9 show the 
same figures with numerical basins of attraction superimposed for four different At by 
the R-K 4 and implicit Euler methods, respectively. Note how little information 
Figures 6.6 and 6.7 can provide as compared to Figures 6.8 and 6.9. Note also how 
rapidly the size of the basin (red) decreases as At increases for the R-K 4 method. This 
phenomenon can relate to practical computations where only a fraction of the 
allowable linearized stability limit of At is safe to use if the initial data is not known. For 
At = 1.75 and 2, spurious limit cycles of higher order period exist. (The multiple white 
circles with only one distinct basin of attraction). In this case, the red regions represent 
the basins of the spurious numerical solutions. 

Figure 6.9 illustrates the situation where unconditionally stable LMM schemes can 
converge to a wrong solution if one picks the initial data inside the green region which 
are valid physical initial data for the ODE. Thus even though LMM preserved the same 
number of fixed points as the underlying ODE, these fixed points can change type 
and stability. This phenomenon is related to the “non-robustness” of implicit methods 
sometimes experienced in CFD computations. In this type of computation where the 
initial data are not known, the highest probability of avoiding spurious asymptotes is 
achieved when a fraction of the allowable linearized stability limit of At is employed. 
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Figure 6.3 (See Color Plate III at the back of this issue.) 


6.3 Numerical Results for the Damped Pendulum Equation 

Selected results for the studied numerical methods for r. = 1 and t: = 1.5 are shown in 
higutes 6. 10 6.12, Here, for each At value the different colors represent different 
numerical basins of attraction of the respective asymptotes. Observe the striking 
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Figure 6.4 (See Color Plate IV at the back of this issue.) 


difference in behavior between the explicit and implicit methods. The shapes and sizes 
of the numerical basins of attraction by the implicit Euler method for At = 0. 1 shown in 
Figure 6. 1 2 appear to be similar to the exact basins of attraction of the DE (4.2). From 
the different colors of the basins in Figure 6.10 one can readily identify that spurious 
higher than period one and spurious limit cycles exist for the different At values by the 
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Figure 6.5 (See Color Plate V at the back of this issue.) 
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Figure 6.6 Phase Trajectories Dissipative Complex Equation, t: = 1, R-K 4. 


explicit Euler method. For At = 1.4, the explicit Euler produces spurious period two 
fixed points. Figure 6.11 shows the existence of spurious fixed points below the 
linearized stability limit by the modified Euler method and spurious fixed points of 
period 4 (the four white dots one each basin) above the linearized stability limit by the 
R-K 4 method. Figure 6.12 shows the evolution (birth and death) of spurious fixed 
points of higher-order period for the implicit Euler method. This figure illustrates 
another situation where unconditionally stable schemes can converge to a wrong 
solution even though these schemes preserved the same number and type of fixed points 
as the underlying ODE. In this case it is the birth of spurious stable and unstable 
asymptotes or even numerical chaos that contributes to the size reduction of the true 
basins of attraction of the ODE. 

6.4 Numerical Results for the Predator-Prey Equation 

Selected results for the two representations of numerical basins of attraction are shown 
in Figures 6.13 6.18. Comparing Figures 6.13, 6.15, 6.16, with Figures 5.4 and 5.7, one 
can again appreciate the added information that the basin of attraction diagrams can 
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Figure 6.7 Phase Trajectories Dissipative Complex Equation, e = 1 Linearized Implicit Euler. 


provide. Here for all of these figures (except Fig. 6.15 for the last four At values), the 
green regions represent the numerical basins of attraction for the stable spiral 
(2.1, 1.98) and red regions represent the numerical basins of attraction for the stable 
node (0, 0). 

The numerical basins of attraction in Figure 6. 1 8 with Ar = 0. 1 appear to be the same 
as the exact basins of attraction of the DE (4.3). The numerical basin of attraction by the 
implicit Euler for the fixed point (0,0) with Ar = 0.1 is larger than the corresponding 
exact basin of attraction for the DE (4.3). In this case the numerical basin of attraction 
for the divergent solution (black region) is smaller than the true one. The dramatic 
difference in shapes and sizes of numerical basins of attraction for the different methods 
and solution procedure combinations compared with the exact basin of attraction is 
even more fascinating than for the previous two models. 

Take, for example, one of the most interesting cases, the modified Euler method. 
Figure 6.15 shows how spurious stable fixed points can alter the numerical basins of 
attraction of the stable node and spiral of the ODE (4.3) for time steps that are below the 
linearized stability limit of both of these stable fixed points of the ODE (see Figs. 5.4 



256 


H.C. YEE AND P. K. SWEBY 


Basins of Attraction 
Dissipative Complex Eqn., e=1 
R-K 4 

At = 0.5 At =1.5 



At =1.75 At =2 




i 






f igure 6.8 (See Color Plate VI at the back of this issue.) 

and 6. 1 3). For At = 0.8, the stable node bifurcates into a spurious fixed point. Without 
performing the bifurcation analysis one would not be able to detect this particular 
spurious fixed point, since the value of the spurious one is so close to the exact fixed 
point U E = (0, 0). For At = 0.9524, there is the birth of a spurious limit cycle (the white 
close curve). For At — 1 .2, spurious higher-order periodic solutions exist. Note that for 
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Figure 6.9 (See Color Plate VII at the back of this issue.) 


the first four At values in Figure 6. 1 5, the fixed points and asymptotic values are colored 
black instead of white due to the birth of additional numerical basins of attraction that 
are colored white. 

The implicit methods change the two saddle points into stable or unstable fixed 
points of other types as illustrated in Figures 6. 1 4, 6. 1 7 and 6. 1 8. For the implicit Euler, 
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Figure 6.10 (See Color Plate VIII at the baek of this issue.) 


the two fixed points (2.1, 1.98) and (0,0) are unconditionally stable and the stabilized 
fixed points (1,0) and (3,0) (saddles for the original ODE) are almost unconditionally 
stable except for small At. This is most interesting in the sense that the numerical basins 
of attraction for the stable exact fixed points U E of the model (4.3) by the implicit Euler 
method were permanently altered for At near or larger than 3 as illustrated in 
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Figure 6.11 (See Color Plate IX at the back of this issue.) 
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Kigure 6.12 (Sec Color Plate X at the back of this issue.) 


Figures 6.14,6.1 7. It would be easier to interpret the results in Figure 6.14 if one 
interchanged the yellow and green colors for At ^ 1. Observe how the newly created 
numerical basins of attraction by the stabilized fixed points (1,0) and (3,0) resulted in 
the segmentation of the numerical basins of attraction of the stable node (0,0) and 
stable spiral (2.1, 1.98). Although the trapezoidal method did not turn the two saddle 
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Pigure6.13 (See Color Plate XI at the back of this issue.) 


points (1,0) and (3,0) into stable fixed points of different type, they did turn the two 
saddle points into unstable fixed points of different type. 

The evolution of the numerical basin of attraction as At changes is very traumatic 
these implicit LM Ms. The cause of nonconvergence of these implicit LM Ms may due to 
the fact that their numerical basins of attraction are fragmented. Take for example the 
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Figure 6.14 (See Color Plate XII at the back of this issue.) 
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trapezoidal method (Fig. 6. 1 8) where the scheme becomes effectively unstable for large 
At. The size of the numerical basins of attraction for the stable fixed points U E shrink to 
almost nonexistence. This phenomenon might be one of the contributing factor to the 
unpopularity of the trapezoidal method in CFD. The basins are so fragmented and 
small for large At that they are beyond the accuracy of the CM2 to resolve and no 
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Figure 6.15 (See Color Plate XIII at the back of this issue.) 
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Figure 6.15 (Continued) (See Color Plate XI If at the back of this issue.) 

further attempt was made. A better approach in computing these types of basins is to 
use interval arithmetic or the enclosure type method (Adams, 1990). 

6.5 Numerical Results for the Perturbed Hamiltonian Equation 

Selected results for the two representations of numerical basins of attraction of the 
various numerical methods for j; = 0.1 are shown in Figures 6.19 6.25. Our studies 
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Figure 6.16 (See Color Plate XI V at the back of this issue.) 


indicate that all of the studied Runge-Kutta methods exhibit spurious limit cycles and 
other spurious periodic solutions. For the Kutta and Heun methods, stable spurious 
asymptotes can occur below the linearized stability limit of the scheme. The implicit 
methods also exhibit spurious asymptotes. I n particular, unstable spurious asymptotes 
were produced below the linearized stability limit by all of the studied schemes. 
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Although this example consists of an artificially small number of grid points, it can 
shed some light on the interplay between initial data, spurious stable and unstable 
asymptotes, basins of attraction and the time-dependent approach to the asymptotic 
numerical solutions. A solid understanding of this concept at the fundamental level can 
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Figure 6.17 (See Color Plate XV at the back of this issue.) 
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Figure 6.17 ( Continued ) (See Color Plate XV at the back of this issue.) 
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help to determine the reliability of the time-dependent approach to obtaining steady- 
state numerical solutions. 

In all of Figures 6.19 6.25, red regions represent the numerical basins of attraction 
for the stable spiral (1/3, 1/3) when At is below the linearized stability of the scheme. 
When At is above the linearized stability, some of the red regions represent the 
numerical basin of attraction of the stable spurious asymptotes. The numerical basins 
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Figure 6.18 (See Color Plate XVI at the back of this issue.) 
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Figure 6.18 ( Continued ) (See Color Plate XVI at the back of this issue.) 
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Figure 6.19 (See Color Plate XVII at the back of this issue.) 


of attraction in Figures 6.21 with At = 0.1 appear to be the same as the exact basins of 
attraction. Note also that the possibility of the numerical basin of attraction being 
larger than the exact one does not always occur when the time step is the smallest. The 
numerical basin of attraction for (1/3, 1/3) is larger than the corresponding exact basin 
of attraction for Af = t by the improved Euler and Kutta methods and for A/ = 0. 1 by 
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Figure 6.20 (See Color Plate XVIII at the back of this issue.) 
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Figure 6.21 (See Color Plate XIX at the back of this issue.) 

the implicit Euler and trapezoidal methods. See Figures 6.21 6.24. The following 
discusses results for the improved Euler, the Kutta, the implicit Euler and the 
trapezoidal methods. 

Improved Euler Method : This example illustrates the existence of spurious limit cycles 
and its effect on the numerical basins of attraction for the exact steady stale. Figure 6.2 1 
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Figure 6.22 (See Color Plate XX at the back of this issue.) 


shows the basins of attraction of the improved Euler method for 4 different At = 0. 1 , 1 , 
2.25,2.35 with /; — 0.1. By a bifurcation computation shown in Figure 6.19, we found 
that the first two time steps are below the linearized stability limit around the exact 
stable steady state (1/3 13), and the last two time steps are above the limit. 
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Above the linearized stability limit spurious limit cycles and higher dimensional 
periodic solutions were observed. Further increasing Af resulted in numerical chaos- 
type phenomena and eventually divergence (with additional increase in Af). For 
At = 2.25 and 2.35, the red or multicolor regions are the basins of the spurious limit 
cycle (the irregular white closed curve shown on Fig. 6.21) or other type of spurious 



Figure 6.23 {See Color Plate XXI at the back of this issue.) 





DYNAMICAL STUDY OK SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 275 



Figure 6.23 [Continued) (See Color Plate XXI at the back of this issue.) 
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Figure 6.24 (Sec Color Plate XXII at the back ol this issue.) 


asymptote (white dots for Kig. 6.21). For these two time steps the numerical basins for 
theexact steady state (1/3, 1/3) by the improved Euler method disappeared. However, if 
the initial data are in the red or multicolor region, one gets spurious solution instead of 
what the linearized stability predicts divergent solution. 
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Figure 6.25 (See Color Plate XXIII at the back of this issue.) 
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Kutta Method : To give an example of the existence of spurious stable asymptotes below 
the linearized stability limit of the scheme, as well as the existence of spurious limit 
cycles above the linearized stability limit. Figure 6.22 shows the basins of attraction 
of the Kutta method for 6 different fixed time steps At = 0.1, 1, 1.826, 1.85,2.75 and 
2.785 (the first four below the linearized stability of the scheme) with £ = 0.1. For 
At = 1.826, the numerical basin for (1/3, 1/3) has become fractal like with the birth of 
fragmented, isolated new basins of attraction due to the presence of spurious periodic 
solutions (the three white complicated closed curves with the associated purple, green 
and blue basins shown in Fig. 6.22). The last two time steps in Figure 6.22 show the 
diappearance of the numerical basin of attraction for the exact steady state with the 
birth of basins for the spurious limit cycle. 

Implicit Euler Method: This is yet another interesting illustration of the use of an 
unconditionally stable implicit method where in practical computations, when the 
initial data are not known, the scheme has a higher chance of obtaining a physically 
correct solution if one uses a At restriction slightly higher than that for the stability limit 
of standard explicit methods (but with larger numerical basins of attraction than the 
explicit method counterparts). Figures 6.20 and 6.23 show the two representations of 
numerical basins of attraction using the implicit Euler method. These figures show the 
generation of stable spurious asymptotes for At ^ 1. As At increases further, the size of 
the same numerical basin decreases and becomes fractal like, and new numerical basins 
are generated. The behavior is similar to the predator-prey model (4.3) in a sense that 
the numerical basin of attraction for (1/3, 1/3) was permanently altered for At near or 
larger than 10. Observe the fragmentation of the numerical basin of attraction for 
(1/3, 1/3) by the basins of the spurious asymptotes. 

Trapezoidal Method: Figures 6.20 and 6.24 show the two representations of numerical 
basins of attraction using the trapezoidal method. As in the implicit Euler case, this 
scheme has a higher probability of obtaining a physically correct solution if one uses a 
At similar to that of standard explicit methods (but with larger numerical basins of 
attraction than the explicit method counterparts). In a manner similar to the implicit 
Euler, the numerical basins of attraction for (1/3, 1/3) are much larger than the 
corresponding exact basin of attraction for At ^ 2. Their sizes are bigger than the ones 
generated by the implicit Euler method with the same At values. The scheme becomes 
effectively unstable due to the fragmentation of the numerical basins of attraction. 
Again due to the high cost of double precision computations, no further attempts were 
made for Ar large. The computation of these basins requires an interval arithmetic or 
the enclosure-type (Adams, 1990) of mathematical operation before a more precise 
behavior can be revealed. 

Straight Newton vs. Other Studied Methods: Figure 6.25 shows the basin of attraction 
using Newton method in solving the steady part of the ODE (hereafter referred to as 
straight Newton) compared with the implicit Euler at Ar = 1 . One can see that straight 
Newton method has a smaller attracting basin for the stable spiral (1/3, 1/3) than the 
implicit Euler method for small Ar. In fact its basin is the same as the implicit Euler 
using larger Ar. Figure 6.25 illustrates the situation where quadratic convergence by the 
Newton method can be achieved only if the initial data are in the red regions. 
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Figures 6.23 and 6.25 also illustrates the fact that using very large At by the (linearized) 
implicit Euler method has the same chance of obtaining the correct steady state as the 
Newton method if the initial data are not known. Comparison of Newton method with 
other iteration procedures for the implicit Euler and trapezoidal methods are reported 
in our companion paper (Yee and Sweby, 1993a). 

Combining the current result with Yee and Sweby (1993a), we can conclude that 
contrary to popular belief, the initial data using the straight Newton method may not 
have to be close to the exact solution for convergence. Straight Newton also exhibits 
stable and unstable spurious asymptotes. Initial data can be reasonably removed from 
the asymptotic values and still be in the basin of attraction. However, the basins can be 
fragmented even though the corresponding exact basins of attraction are single closed 
domains. The cause of nonconvergence may just as readily be due to the fact that its 
numerical basins of attraction are fragmented. 

6.6 Global Asymptotic Behavior of Iterative Implicit Schemes 

The global asymptotic nonlinear behavior of some standard iterative procedures in 
solving nonlinear systems of algebraic equations arising from four implicit linear 
multistep methods (LM Ms) in discretizing models (4. 1 ), (4.3) and (4.4) is analyzed in our 
companion paper (Yee and Sweby, 1993a). The implicit LMMs include implicit Euler, 
trapezoidal, mid-point implicit and three-point backward differentiation methods. The 
iterative procedures include simple iteration and full and modified Newton iterations. 
The results are compared with standard Runge-Kutta explicit methods, a non-iterative 
implicit procedure, and straight Newton method. Here we give a summary of Yee and 
Sweby (1993a) so that the reader may get a bigger picture of implicit methods other 
than the ones studied in this paper. 

Studies in Yee and Sweby (1993a) showed that all of the four implicit LMMs 
exhibit a drastic distortion but less shrinkage of the basin of attraction of the true 
solution than standard explicit methods studied in this paper. In some cases with 
smaller A;, the implicit LMMs exhibit enlargement of the basins of attraction of 
the true solution. Overall, the numerical basins of attraction of a non-iterative im- 
plicit procedure mimic more closely the basins of attraction of the continuum than 
the studied iterative implicit procedures for the four implicit LMMs. In general the 
numerical basins of attraction bear no resemblance to the exact basins of attraction. 
The size can increase or decrease depending on the time step. Also the possible existence 
of the largest numerical basin of attraction that is larger than the exact one does 
not occur when the time step is the smallest. The dynamics of numerics of the implicit 
methods differ significantly from each other, and the different methods of solving 
the resulting non-linear algebraic equations are very different from each other 
since different numerical methods and solution procedures result in entirely different 
nonlinear discrete maps. Although unconditionally stable implicit methods allow 
a theoretically large time step At, the numerical basins of attraction (allowable initial 
data) for large At some-times are so fragmented and/or so small that the safe (or 
practical) choice of Af is slightly larger or comparable to the stability limit of standard 
explicit methods (but with larger numerical basins of attraction than the explicit 
method counterparts). I n general, if one uses a At that is a fraction of the stability limit. 
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one has a higher chance of convergence to the correct asymptote than the standard 
explicit methods. 

Studies in Yee and Sweby (1993a) also showed that the variable time step control 
method can occasionally stabilize unstable fixed points, depending on the initial data, 
starting time step and the iterative tolerance value. One shortcoming is that the size of 
At needed to avoid spurious dynamics is impractical to use, especially for the explicit 
method. 


7. CONCLUDING REMARKS 

The global asymptotic nonlinear behavior and bifurcation phenomena for the explicit 
Euler method, five different multistage Runge-Kutta methods (modified Euler, improved 
Euler, Heun, Kutta and 4th-order methods), two and three-step predictor-corrector 
methods, Adams-Bashforth method, and implicit Euler and trapezoidal method with 
linearization are compared for different model nonlinear ODEs. The five multistage 
Runge-Kutta methods and the predictor-corrector methods are nonlinear in the 
discretized parameter space At and all LMMs are linear in At. With the aid of the 
CM-2, the complex behavior and sometimes fractal like structure of the associated 
numerical basins of attraction of these time discretizations are compared and revealed 
for the first time. 

The numerical results indicate that with sufficiently small Af and initial data close to 
the steady state (usually not known for the time-marching method), one can have the 
highest chance of convergence to the correct asymptote. In general, the initial data can 
be far removed from the exact steady state by the studied implicit methods provided 
that a fraction of the allowable time step restriction is used. Our study also indicates 
that bifurcation to a period two or lower order period soluton is readily detectable in 
numerical calculations. However, bifurcation to a limit cycle will not be so obvious 
(without a phase portrait representation), especially in the vicinity of the bifurcation 
point. Indeed the phenomenon of an artificial time iteration to steady-state of a large 
system formed by spatial discretization which nears convergence before the residuals 
"‘plateaus out", could actually be the result of a stable spurious limit cycle around the 
Hopf bifurcation point. In addition, the bifurcation of spirals to limit cycles might account 
in part for the phenomenon of near (but lack of ) convergence in large stiff systems. 

For a given initial data and two finite but different At's that are below the linearized 
stability limit of the scheme, their numerical solutions might converge to two different 
solutions even if no spurious stable steady-state numerical solution is introduced by the 
scheme and the initial data are physically relevant. The source of the behavior is due to 
the existence of unstable spurious asymptotes or stable asymptotes other than steady 
states which have the same detrimental (in terms of robustness) effect. However, in the 
case of occurrence of stable spurious steady states, they can be mistaken for the true 
steady state in practical computations. In other words depending on the initial data, for 
a given At below the linearized stability limit, the numerical solution can (a) converge 
to the correct steady state, (b) converge to a different steady state, (c) converge to 
a spurious periodic solution, (d) yield spurious asymptotes other than (a)-(c), or (e) 
diverge, even though the initial data are physically relevant. 
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Another important finding is that unlike the scalar first-order autonomous ODE 
discussed in part I (Yeeef al., 1991), the fixed points can change types as the time step is 
varied even for two-time-level unconditionally stable implicit LMMs. An unstable 
fixed point can become a stable fixed point and can e.g., change from a saddle to a stable 
or unstable node (for a fixed system parameter /:). Since these implicit methods can 
introduce spurious asymptotes as well, thus even though LMMs preserve the same 
number but not the same types of fixed points as the underlying DEs, the numerical 
basins of attraction of LM Ms (explicit or implicit) do not always coincide with the exact 
basins of attraction of the underlying DEs. One major consequence of this behavior is 
that the flow pattern can change type as the discretized parameter is varied. Another 
consequence of these phenomena is the fragmentation of the numerical basin of attrac- 
tion. In general, unconditionally stable implicit LMMs exhibit less shrinkage of the 
basin of attraction of the true solution than standard explicit methods. Another 
interesting result is that contrary to popular belief, the initial data using the straight 
Newton method may not have to be close to the exact steady state for convergence. 
However, we believe that one cause of nonconvergence in straight Newton or implicit 
LMMs with large time step may due to the fact that the numerical basins of attraction 
are fragmented. 

In conclusion, the present results can explain some of the roots of why one cannot 
achieve the theoretical linearized stability limit of the typical implicit Euler and 
trapezoidal time discretization in practice when solving strongly nonlinear DEs, e.g. in 
CFD. The results can also shed some light in bridging some of the gaps between 
theoretical convergence criterion (At ^0, as n -* oc) and practical scientific computa- 
tion (finite At as n -» oc). 
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Figure 6.1 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.2 (See H.C. Yce and P. K. Sweby.) 
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Figure 6.3 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.4 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.8 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.10 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.12 (See H.C. Yee and P. K. Sweby.) 
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Color Plate XI Figure 6.13 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.15 (Continued) (See H.C. Yee and P. K. Sweby.) 



Basins of Attraction 
Predator-Prey Eqn. 


R-K 4 


At =0.4 


At = 0.8 


At = 1 


At = 1.77 


Color Plate XIV Figure 6.16 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.17 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.19 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.23 (See H.C. Yee and P. K. Sweby.) 
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Figure 6.24 (See H.C. Yee and P. K. Sweby.) 
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Figure 10 Computed surface pressure contours with store located in captive position 
and moving through 0.6, 1 .0, 1.6, and 2.0 store diameters below the pylon. 
(See A. Arabshahi and D. L. Whitfield) 


